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"Method for determining the three-dimensional surface of an object". 

* * * * 

DESCRIPTION BACKGROUND OF THE INVENTION 
Field of the Invention 
5 The present invention refers to a method for determining the three- 

dimensional surface of an object and a_relevant computer program. 

Description of Related Art 

For describing and analysing the geometry of objects that have to be 

10 produced, modelling methods that use computers are known. In particular, 

when a three-dimensional surface enclosed by a set of points acquired from 
a real object has to be obtained, a technique called level-set that describes an 
enclosed surface as a level surface of a volumetric function (or implicit 
function) is normally used. 

1 5 The level-set technique provides for the time evolution of the 

volumetric function, according to a suitable equation atof the partial 
derivatives, typically belonging to the Hamilton - Jacobi family of equations. 

The volumetric function describes the distance of each point from the 
enclosed surface and this distance will have positive or negative sign 

20 according to the point considered whether it is inside or outside the surface. 

Thus the set of points that possess a null distance from the surface S 
represent the surface itself. The technique provides for the evolution of the 
enclosed surface, defined by the Hamilton - Jacobi equation so that it 
surrounds and adheres to the cloud of 3D points that describes the object. 

25 This technique is used as it is capable of adapting to the various types of 

objects and to the most diverse point acquisition techniques. This technique 
in fact, does not require any additional information about the topology of the 
points acquired and is insensitive to the technique used: the density of the 
points can vary in the different zones of the object considered; no 

30 information is required about the point acquisition sequence or about their 



relative spatial positions; and the number of objects enclosed separated from 
each other is not relevant for the algorithm. 

As this technique requires the updating of a volumetric function at each 
step of evolution of the front, the calculation cost is usually very high, so 
much so that it makes the method of very little practical utility. Methods are 
known that permit the reduction of the calculation cost of this technique, for 
example limiting the updating of the volumetric function in an enclosed 
zone surrounding the front in evolution. Despite this the calculation cost 
remains prohibitive. 

T e chniques in alt e rnativo Alternative techniques to the level-set 
technique are also known. For example, methods such as those of the 
Delaunay triangulation are known, as well as the NURBS (Non Uniform 
Rational B-Splines) and the HRBF (Hierarchical Radial Basis Function). 
Each of these methods has limitations and impes eimposes a series of 
constrictions restrictions on the set of the points acquired, such as for 
example the calculation complexity depends heavily on the number of 
points, or they need to have a spatial distribution of the points acquired as 
uniform as possible. 

BRIEF SUMMARY OF THE INVENTION 
In view of the state of the technique described, an object of the present 
invention is to provide for a rapid method for determining the three- 
dimensional surface of an object and-that does not have the inconveniences 
of the known art. 

In a first aspect^ the present invention relates to a method for 
determining the three-dimensional surface of an object comprising the 
phases of: defining the coordinates of a plurality of points of said object; 
defining a three-dimensional matrix of cells that contains said object to 
which a value can be associated; determining the distance between each 
centre of said cells of said three-dimensional matrix of cells and the closest 
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point of said plurality of points of said object; setting the value of several 
cells of said three-dimensional matrix of cells at a first preset value; 
determining the value that each cell of said three-dimensional matrix of cells 
assumes, with the exception of said several cells, by means of the following 
5 formula 

F(x l ,0-v l +w-£F(x / ,0-v, 
F(x,,r + 1) = ± 

j 

where 

5t, represents the coordinates of the centre of the i_th cell, 
F(x jy t) represents the value of the i_th cell at time t, 
v i represents said distance, 
10 w represents a second preset value, and 

j indicates a neighbourhood of cells of the i_th cell; 
determining the module sum of the variations of the value of each cell 
between the time t and the time t+l\ and repeating said phase of determining 
the value that each cell of said three-dimensional matrix of cells assumes if 
1 5 said sum is higher than a third preset value. 

In its second aspect^ the present invention relates to a computer program 
comprising a program code that carries out all the phases of the method for 
determining the three-dimensional surface of an object when said program is 
carried out on said computer. 
20 In its third aspects the present invention relates to a computer program 

memorised on a suppor t stored in a memory that can be used by saida 
computer to control the execution of all the phases of the method for 
determining the three-dimensional surface of an object. 

Thanks to ln the present invention, by substitutin g replacing the 
25 differential equation atof the partial derivatives normally used in literature 

(Hamilton-Jacobi) with one that is very well suited to describing 
discontinuous functions such as those that describe the interface between 



two immiscible fluids, that is one of the Navier-Stokes equations, derived 
from the fluid dynamics, the calculation cost is considerably reduced, and a 
much quicker convergence can be obtained. The surface can be interpreted 
as the front between two viscous fluids that moves until it completely 
5 encloses the object. In particular, acting on the viscosity value of the fluids 

used it is also possible to control the level of roughness permitting smooth 
surfaces to be obtained even if the points presented deviations from their 
correct spatial position. In addition, the method proposed permits the 
introduction of elements and strategies that enable an end control of the 

10 evolution of the front, based on simple principles which are interpretable 

phvsicallv of simple physics intorprotabilitv . For example, an artificial 
turbulence can be introduced in the fluids that facilitates penetration inside 
particularly narrow cavities. Thus it is possible to provide a correct 
determination even of configurations of points that traditionally had been 

15 negated. Thanks to its flexibility, the Navier-Stokes equation can be used for 

determining more volumetric functions that describe different parts of the 
object considered. Objects that are particularly complex or that have high 
resolutions can thus be represented by calculating separately the evolution of 
the volumetric function on the different parts that make them up. The 

20 siftg teindividual results obtained are then recomposed to obtain the overall 

volumetric function using a simple mathematic operator. 

BRIEF DESCIRPTION OF THE DRAWINGS 
The characteristics and advantages of the present invention will appear 
25 evident from the following detailed description of an embodiment thereof, 

illustrated as non-limiting example in the enclosed drawing, in which Figure 
1 shows a flow diagram of the method for determining the three-dimensional 
surface of an object. 
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DETAILED DESCRIPTION OF THE INVENTION 
The present invention is based on the Navier - Stokes differential 
equation at the partial derivatives for the conservation of mass in which a 
redefinition of the velocity vector is used according to a criterion described 
5 hereunder. It describes the motion of two immiscible fluids that converge, 

one from inside and the other from outside, towards the cloud of points that 
describes the object considered. At the end of the evolution, when both the 
fluids adhere to the points of the object, the surface of separation between 
the two coincides with the surface of the object itself. 

10 The volumetric function describes the contents of the two fluids in each 

elementary cell or voxel, where voxel means the elementary element of the 
volumetric grid, such as how the pixel is used for a two-dimensional grid. 
To distinguish the two fluids we have assigned opposite sign mass to them. 
Therefore a cell in which the volumetric function is worth -1 contains 

15 exclusively the internal fluid, while a value of +1 indicates that only external 

fluid is present. 

If between two adjacent cells the volumetric function changes sign it 
implies that the surface of separation passes between them and thus also the 
surface of the object. The exact point of passage is determined by means of 
20 interpolation operations. 

The exact position of the contact surface between the two fluids is thus 
identified by the position of the level-set zero (that is the set of points with 
null distance from the surface). 

The law for the conservation of mass is one of the fundamental 
25 principles of elassi eclassical physics and it is independent from the nature of 

the fluids and from the forces that act on it; it expresses the principle 
according to which in a fluid dynamic system mass cannot disappear or be 
created unless sinks or sources are present. 

The value F of the volumetric function in each cell indicates the 
30 quantity of fluid inside it and this value is modified by the forces acting on 



the fluid, by its velocity, by the diffusion coming from the adjacent cells and 
from the sources or sinks. 

The evolution of a fluid inside a system can be described by means the 
flow vector G that expresses the quantity of fluid that crosses a surface in a 
given interval of time. 

In its classic formulation the law of conservation for mass affirms that 
the time variation by time unit of the quantity of matter F inside a volume Q, 




is equal to the net contribution of the flow that crosses the external 
surface S of Q, to the contribution Qy of sources (or sinks) of volum e inside 
the volume -it and to the contribution of sources O s oriented oricntatod, along 
its surface S. The contributions due to the different sources can thus be 
written as: 



The general formulation for the conservation of mass equation can thus 



be written as: 

The negative sign in front of the flow integral is due to the fact that the 
area dS is orientated outwards from Q. 

Using the Gauss theorem the previous equation can be rewritten using 
exclusively volume terms: 



The flow density vector G can be decomposed in a diffusive 
component G D and in a convective component G c , the convective part of 
the flow vector describes the transport phenomenon due to the external 
forces and is defined as the product between the velocity v of the flow and 



\Q v dn+\Q s dS 



n s 
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the quantity of matter transported F: therefore it results G c =vF . 

The diffusive component instead derives from the remixing, present also 
in resting fluids, due to the-thermal agitation, and is usually proportional, 
according to the law of Fick, to the gradient of F 

G D = f7F where /is the diffusivity constant. The time evolution of the 
volumetric function F due to the diffusivity is regulated (as can be verified 
forcing the conservation of mass) by the differential heat equation at th e 
partial derivatives of the heat and thus its contribution results are mainly 
connected to regularizing and a flattening of the function F. The effect of the 
diffusive contribution is thus taken into consideration applying a Gaussian 
filtering to the entire volumetric function. 

The intensity and the orientation of the velocity v of the fluid is 
determined in each cell of the matrix. This velocity has been defined 
orientating the flow towards the point of the 3D cloud (that describes the 
object) closest to the cell, and the module modulus of this vector has been 
selected directly proportional to the distance: in this manner both the fluids 
progressively slow down coming closer to the cloud of points and converge 
at the final surface of contact. The equation that has thus been used to 
describe the evolution of the two fluids is the following: 



The local variation of the quantity of fluid contained in the i_th cell 
F(x i9 t) between two successive instants therefore results determined by the 



Gaussian function whose standard deviation a describes the degree of 
diffusion of the fluid between the different cells. 

The velocity v(x) has been defined for each cell according proport 




convective contribution deriving from each cell F(x,/)|v(x)| weighed with a 
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to the distance between the centre of the cell and the closest point belonging 
to the 3D cloud. In particular the module modulus of this velocity equals 
|jc,. ~p\ a where x i is the centre of the i_th cell, while p is the closest point 
of the cloud of points (the point at least distance from the centre of the i _th 
5 cell), a is instead a parameter that regulates the course of the velocity 

function. Experimentally it has been determined that acceptable values for a 
are between 4r£1.5 and 2r#2.5. The greater this value is, the slower the 
velocity will be in proximity efto the points while it will be rapid at a great 
distance. Even though this characteristic is very important for a rapid 

10 evolution of the values of the matrix it can however lead to a slow final 

convergence. A preferred value of the value of a equals 2.0. 

Another property of the fluids used in the model is the viscosity, whose 
value determines the standard deviation a inside the equation. This property 
describes the viscos e viscous friction that acts inside the fluid limiting its 

1 5 mobility. Considering this property we prevent the fluids that cross from 

crossing the reconstructed surface from and flowing between the interstices 
of the various points. This property is particularly important as through it 
special arrangement of points can be correctly met such as fine blades or 
particularly pointy objects. Thanks to the viscosity, the fluids do not cross 

20 the surface which remains adherent to the blade or to the point. It is also 

possible to reduce the curving of the entire surface obtaining particularly 
smooth surfaces and reducing the surface roughness that can also have 
originated from the error in measurement. It must be remembered however 
that high values can how e ver limit the capacity of the algorithm to represent 

25 objects with particularly jagged surfaces and to enable the flow of the liquid 

inside concave zones. Because of the above-mentionedi the viscosity is a 
parameter whose choice is determined according to the typology of the 
object considered, by the imprecision of the acquisition system used A and by 
the degree of definition required for the adhesion of the final surface to the 

30 cloud of points. 



Another aspect that has been used fey]n the physical modelling of the 
fluids is the-turbulence, which treats the presence of whirlpools that describe 
local variations of the velocity function. 

The presence of the-turbulence can locally develop high pressures and 
can generate instability in the flow but at the same time can favour the 
penetration of the fluid inside very small openings and prevent stagnation in 
zones with modest velocities. 

The presence of vortexes wasis introduced by modulating the field of 
velocity with small oscillations in the different directions around the initial 
configuration. Because of the close correlation between the velocity vector 
in each cell and the vector that points towards the closest point of the 3D 
cloud, the local oscillations of the velocity field can also be interpreted as an 
elastic deformation of the cloud of points in the various directions which can 
facilitate the entrance of the fluid inside the restricted zones. 

The modulation of the velocity function is applied cyclically along the 
three spatial axes by multiplying by a factor the terms v of the velocities of 
the cells that lie orientated in parallel to these axes in relation to the i_th cell 
under examination and its intensity can be controlled externally, in particular 
at each iteration the velocity is amplified in a given direction and/or the 
velocity in the other two directions is reduced by an equal quantity. 
Obviously the g reater the expansion coefficients used the greater will be the 
turbulence. Values that are too high can however determine an excessive 
remixing alse-between the two fluids on the contact surface causing the 
leakage of the internal fluid or an excessive penetration of the external fluid 
in the object considered. 

The previous equation that was used to describe the evolution of the two 
fluids has been discretized to be able to carry out the operations more easily 
by means of a computer and is described by the following equation. 



10 



j 

where 

3c, represents the spatial coordinates (three-dimensional array) of the 
i_th cell, 

F(x i9 t) represents the value of the i_th cell at time /, 

v, represents the velocity (or the distance as previously defined), 

w represents a preset value of viscosity, and 

j is a preset value that indicates a neighbourhood of cells of the i_th 

cell. 

In particular, the index summation j considers the contribution of the 
closest cells to the point under examination. Preferably, this interaction was 
chosen to be limited exclusively to the first 6 cells near each cell under 
examination, that is that share with it an entire face of the cubical cell (the 
cells at distance 1 from the cell under examination). 

We now refer to Figure 1 that shows the flow diagram of the method for 
determining the three-dimensional surface of an object. This flow diagram 
describes how a computer program carries out the calculation for 
determining the three-dimensional surface of an object. 

Phase 1 - Data definition. 

In this phase the coordinates of a plurality of points of the object that is 
needed to be represented are defined. The set of the points is described as a 
set of vectors with 3 components, containing the spatial coordinates of each 
point. No additional information is requested in relation to the points. 

Phase 2 — Definition of cell matrix. 

In this phase the cell matrix (or three-dimensional grid) that defines the 
evolution dominium domain of the value of F is described. The number of 
cells is chosen on the basis of the dimensions of the object that is needed to 
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be represented. The number of cells that will be used is determined on the 
basis of the extensions along the three axes of the object. Thus keeping a 
cubical form for each elementary cell of the matrix, a grid that contains 
e xactly outt he object is defined. It is also possible to define a form of the 
cells that is not cubical, that is a different resolution of the matrix along the 
different axes can be set. This can be very important for objects in which a 
high resolution is necessary only in a dimension such as low relief, building 
faces, etc.. 

It is also possible to define a rotated matrix in relation to the axes 
according to which the object has been acquired and, in particular, align it 
with the main directions of the object itself. Thus it is possible to minimize 
the number of cells used at the same time keeping the same resolution. A 
centre with the appropriate coordinates belongs to each cell, and a value can 
be associated to each cell. 

Phase 3 - Determining the velocity. 

In this phase the field of the velocities that belong to each cell of the 
grid is defined. As we saw previously the velocity field depends closely on 
the distance of each centre of the cells from the closest point of the cloud of 
points. Thus, the field of the distances is a three-dimensional matrix with the 
same dimensions as the grid. For each cell the point of the 3D cloud that is 
the closest (at least distance) is determined and the distance is calculated as 
- p\ a between this point and the centre of the cell. This value is stored in 
memory m e moriz e d in an equivalent cell of the matrix of the distances. 

Take note that it is not necessary to know the direction of the velocities 
but only the modulo modulus . 

It is important to note hewthat, the dependence of the 
calculation calculational complexity of this method on the number of points 
of the 3D cloud is limited exclusively to the determination of the velocity 
field, only once at the beginning of the procedure A and all the evolutive 
phase it is completely independent from the number of sample points used. 
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Phase 4 - Arrangement of the sources. 

Different arrangements of the sources are possible. A first arrangement 
provides for the external sources to be placed exclusively on the edges of the 
matrix (a first portion of the cells, which will thus always be kept at the 
value +1 during the entire evolution) and all the remaining cells instead to be 
filled with external fluid (value equal to -1). A second arrangement instead 
provides also for one or more internal sources to be provided, in this case the 
cells of the edge will always be kept at +1, the value of the internal cells will 
be equal to 0, while the internal sources will always be kept at -1 . This 
method permits an even more rapid evolution but requires initial information 
relative to the arrangement of the additional internal sources that areis not 
easily available. 

Phase 5 - Fixing viscosity and turbulence parameters. 

The viscosity (in the discretized formula represented by w) regulates the 
diffusion of the fluid and thus the weight of the interaction between each 
single cell and its neighbours. The preferred value of w resultsi s between 
OrWJ. and 0^O9. 

The parameters relating to the turbulence indicate how the velocities of 
the fluid, in each cell, are altered at each iteration with the purpose of 
simulating the turbulent effect. In particular it is possible to set for each of 
the 3 axes the coefficient by which the velocity in this direction will be 
multiplied cyclically. The multiplicative factors for the velocity vary 
preferably by a minimum of 0.5 up to a maximum of 1 .5. The use of the 
turbulent regime is optional. 

Phase 6 - Starting up the procedure. 

To accelerate the operations, and optionally, provision is made in this 
field for another three-dimensional matrix with the dimensions of the matrix 
of the cells in which the information regarding the updating or not of the 
different cells is contained. If a cell can be updated its contents will be 
processed by the evolutive algorithm and thus the contents of liquid inside it 
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can vary between two iterations. Cells that cannot be updated are the sources 
(that are always kept at the same value, +1 or -1). During the evolution also 
those cells for which the variation of liquid between two iterations has not 
been relevant (below a certain set threshold) can become non-updateable and 
thus these cells are frozen until the evolutive front approaches them. 
Phase 7 - Updating cells. 

It is the main phase of the method in which the value of the quantity of 
fluid present in each eel! of the three-dimensional matrix of cells is 
determined. 

The use of the Gaussian function inside the integral would provide for 
the analysis at each iteration of the interaction of each cell with all the 
matrix of the cells, which would result in being particularly cumbersome in 
computation computational terms. 

In the case with the discretized formula, it is possible and preferable to 
limit the extension of this interaction exclusively to a neighbourhood of cells 
closer to the i_th cell, that can be for example 6 or 1 8 or 26 (the cell 
considered is at the centre of a matrix 3x3x3), that is those that are 
respectively in face, edge or vertex contact (that is at distance 1, a/2 and V3 
from the cell under examination, in the case of cubic cells) from the i_th 
cell. Preferably the 6 closest cells have been chosen, that is the 6 cells that 
contact the 6 faces of the i_th cell. 

Once the cell has been updated, the variation in relation to the previous 
value is assessed and if this value is lower than a threshold set by the user 
the contents of the cell are frozen and the cell is declared not updateable. 

Phase 8 - Visualization 

It is possible, if it interests, to visualize the level-set zero, that is the 
surface of separation between the two fluids. The algorithm used is that 
commonly known with the name of Marching Cubes (see for example the 
article "Marching cubes: a high resolution 3D surface construction 
algorithm", by William E. Lorensen and Harvey E. Cline, from Computer 
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Graphics, Volume 21, Number 4, July 1987) by means of which it is 
possible to obtain an effective triangulation of the level-set zero and a sub- 
pixel resolution thanks to the linear interpolation carried out on adjacent 
cells. The triangulation thus obtained is represented on video, for example, 
5 by means of an interface program called OpenGL that permits visualization 

in real time and a direct interaction with the user who can interactively 
change at will the view point analysing the evolution of the level-set in real 
time. 

Phase 9 - Determining the flow variations. 
10 At the end of each evolution the overall quantity of fluid that has been 

moved is assessed. 

Phase 10 - Evolution cycle. 

If the overall quantity of fluid that has been moved is slower than a 
preset threshold the evolution is interrupted at phase 11, otherwise it carries 
1 5 out further iteration. 

At phase 1 1 , when the calculation is completed, a three-dimensional 
matrix of points is obtained that defines the values that each point assumes 
the volumetric function. 

It is possible to apply a low-pass filter to this matrix so as to obtain a 
20 smoother reconstructed surface removing possible errors of sampling 

without acting on the viscosity parameters (that could prejudice the correct 
determination of concave zones). The implementation of this filter is the 
same as that to carry out a mono one- and two-dimensional filtering: the 
entire volumetric function is conveyed with a three-dimensional matrix 
25 containing values of the filter. For example it can be made in low-pass filter 

using a matrix 3x3x3 in which the central element has the maximum value 
and the other elements possess a value given by the value of a Gaussian 
function depending on their distance from the centre of the matrix. 

With the present method, in which case it is possible and convenient to 
30 decompose a surface of large dimensions in distinct parts, it is possible to 
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calculate the evolution singularly on each part and recompose then the 
various volumetric functions. In this manner it is also possible to use 
different resolutions for the different parts, then recomposing them on a 
single common matrix. It is preferable, in order to improve the link between 
the different parts, between the different defined zones, that there is 
overlapping, that is that part of the end points belonging to a zone are also 
incorporated in the adjacent zone. Once the matrices of the various zones are 
obtained, they are arranged on a common matrix with resolution equal to the 
maximum resolution used. In the case the contents of a cell of the final grid 
are contemporarily described by two volumetric functions (because they 
overlap) the minimum between all the volumetric functions that describe it 
is chosen as value for this cell. 

The method for determining the three-dimensional surface of an object, 
herein described, can be transcribed in a program code, in a manner well 
known to a technician of the sector, that can be momorized stored on any 
type of memory or support (floppy, CD) and/or can be carried out by a 
computer. 
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"Method for determining the throo dimensional surface of an object". 

* * * * 

ABSTRACT 

Th e pres e nt invention r e f e rs to a method for determining the throo 
5 dimensional surfac e of an object. 

In an embodiment thereof thc A method for determining the three- 
dimensional surface of an object comprises the phases of: defining (1) the 
coordinates of a plurality of points of said object; defining (2) a three- 
dimensional matrix of cells that contains said object to which a value can be 
10 associated; determining (3) the distance between each centre of said cells of 

said three-dimensional matrix of cells and the closest point of said plurality 
of points of said object; setting (4) the value of several cells of said three- 
dimensional matrix of cells at a first preset value; determining (7) the value 
that each cell of said three-dimensional matrix of cells assumes, with the 
1 5 exception of said several cells, by means of the following formula 

F(x„/ + 1) = 



V, + w • 



j 



2>, 



where represents the coordinates of the centre of the i_th cell, 

F(x i9 1) represents the value of the i_th cell at time t, 

v. represents said distance, 

w represents a second preset value, and 
20 j indicates a neighbourhood of cells of the i_th cell; 

determining (9) the sum in module of the variations of the value of 
each cell between the time t and the time t+l\ repeating (10) said phase of 
determining the value that each cell of said three-dimensional matrix of cells 
assumes if said sum is greater than a third preset value. (Fig. 1). 
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